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Abstract. 

The realistic description of correlated electron systems has taken an important 
step forward a few years ago as the combination of density functional methods and 
the dynamical mean-field theory was conceived. This framework allows access to both 
high and low energy physics and is capable of the description of the specific physics 
of strongly correlated materials, like the Mott metal-insulator transition. A very 
important step in the procedure is the interface between the band structure method 
and the dynamical mean-field theory and its impurity solver. We present a general 
interface between a projector augmented wave based density functional code and many- 
body methods based on Wannier functions obtained from a projection on local orbitals. 
The implementation is very flexible and allows for various applications. Quantities like 
the momentum resolved spectral function are accessible. We present applications to 
SrVC>3 and the metal-insulator transition in Ca2- a; Sr :z: Ru04. 



PACS numbers: 71.30. +h, 71.15.Ap, 71.27.+a, 75.25.Dk 



DFT++ method implemented with projector augmented waves 



2 



1. Introduction 

The theoretical description of strongly correlated electron systems is one of the main 
challenges in contemporary condensed matter physics. Since the intricate properties of 
such systems have become accessible experimentally, a unified theory that can account 
for material-specific quantities as well as complicated many-body behaviour is needed. 
Traditionally, the state of the art electronic structure theory that captures material 
specific physics is density functional theory (DFT) (see e.g Ref. [1] for a review). 
This theory can be applied to itinerant systems (e.g. metals) where the kinetic energy 
dominates localization, but it fails for the strongly correlated ones (transition metal 
oxides, Kondo systems, etc.). The incorrect behaviour of the DFT for strongly correlated 
systems can be remedied by augmenting the DFT with a local Hubbard-like interaction. 
The resulting lattice Fermion model with ab initio parameters from DFT is then solved 
by the dynamical mean- field theory (DMFT) (see e.g. Ref. [2] for a review). The 
resulting methodologies are commonly referred to as LDA+DMFT [3, 4] (see Ref. [5] 
for a recent review) or in the static Hartree-Fock limit LDA+U [6, 7]. In such acronyms 
LDA is used synonymously with DFT, since the methods are not specifically limited to 
the LDA functional. In general a method interfacing DFT with an additional many- 
body treatment can be called DFT++ following Ref. [4]. DFT++ combines the ab initio 
capabilities of DFT to include material specific effects (crystal structures, etc.) with 
the capabilities of model Hamiltonians to account for the physics of strongly correlated 
electrons (Mott transitions, Kondo physics, etc.). From a practical point of view, all 
DFT+- 1- methods base on a Hubbard type interaction being added to the Kohn-Sham 
Hamiltonian [8], which serves as bare, "non-interacting" starting point. 
An important point at the interface between DFT and many-body methods is a 
suitable choice for a basis of the one-electron Kohn-Sham states. The local Hubbard- 
like interaction in the correlated subspace of orbitals (usually d or / orbitals) is 
best represented in a set of localized orbitals. Thus, earliest implementations of the 
DFT+DMFT used the linear-muffin-tin-orbital (LMTO) basis and represented the 
correlated subspace with the subset of LMTO's with the specific character. This choice 
is certainly sensible, however, it might not be optimal and other basis sets have been 
investigated. A basis set, that has been heavily used in recent years, is the basis of 
Wannier functions. These have been utilized in different flavours in the context of 
DFT+DMFT: Anisimov et al [9], Amadon et al [10], Haule et al [11] and Aichhorn et 
al [12] used different schemes for projections onto a subset of Bloch functions, Pavarini 
et al [13] used Nth order muffin-tin orbitals (NMTOs) and Lechermann et al [14] used 
maximally localized Wannier functions (MLWFs) [15]. 

In this paper we lay out methodological details concerning the use of projections of 
Bloch states onto local orbitals, as introduced in Refs. [10, 16] in the framework of the 
projector augmented wave method [17]. We will refer to the method in what follows as 
PLO (projected local orbitals), as in Ref. [10]. The purpose of this paper is to present 
an implementation of a general DFT++ method in the PAW basis set using the VASP 
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Table 1. Intersite hopping integrals ty%* for SrVC>3 from PLO compared with 
maximally localized Wannier functions in different basis sets and Nth order muffin- 
tin orbitals. Energies in meV. (f from Ref. [14]) 
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[18, 19] code and to explore different methods to construct the underlying Wannier 
functions. We compare the present implementation with maximally localized Wannier 
Functions produced in different basis sets and NMTO calculations. 

The paper is organized as follows: In Section 2 we present a general outline of 
the approach, with special attention to specific details concerning the PAW method. 
The methodological details will be illustrated using SrV0 3 as a test system. As a more 
challenging application of the DFT++ approach a brief discussion of the metal-insulator 
transition in Ca2- :r Sr. r Ru04 is included. 

2. DFT++ in the PAW framework 

2.1. General Presentation 

The ab initio treatment of correlated electron systems requires the calculation of Green 
functions and hybridization functions in terms of local orbitals. This is readily achieved 
when using a basis set, which is localized in real space, such as linear muffin-tin orbitals 
or Gaussian basis sets. However, many implementations of the density functional theory 
use a delocalized plane wave basis set. This has the advantages, that the basis set is 
simple, universal and its convergence is controlled in principle by a single parameter, the 
energy cutoff. The projector augmented wave method (PAW), being a representative 
of plane- wave based methods, is a fast and accurate way of implementing DFT [17]. It 
formally is a so-called all electron method meaning it takes into account all electrons 
in the problem. At the heart of the method is the following theorem, stating that the 
wave function |\&) can be decomposed exactly as follows 

l*> =T|*>= (l + $>)|¥> 

i 

= i*>+E(m w 

i 

Here T is the operator denoting the transformation, which is different from the identity 
only inside an augmentation region, is the so called pseudo wave function and the 
sum runs over all augmentation channels i. As in the construction of pseudopotentials, 
physical partial waves \<f>) are solutions of the Schrodinger equation of isolated atoms, 
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while the corresponding auxiliary (pseudo) functions \<p) are chosen to match |0) outside 
the augmentation spheres, being smooth inside and continuously differentiable in all 
space. The projectors \pi) are finally defined by 

(Pi\4>j) = Sij, 

where the tilde as usually in the PAW formalism discriminates pseudo from physical 
quantities. The Kohn-Sham equations in PAW representation are obtained by applying 
the variational principle to the total energy functional with respect to the auxiliary wave 
functions: Since the transformation operator T does not depend on the electron density, 
the Kohn-Sham equations transform according to 

'PH KS T\V k )=e^T\^) 

Here, Hks is the Kohn-Sham Hamiltonian, so that above equation is a Schrodinger type 
equation, but with the overlap operator occurring on the right hand side. To solve the 
equation the auxiliary wave functions are expanded in terms of plane waves: 

* k (r) = (r|tf k ) = ^c kjG exp(i(k+G)r). 

G 

Following Ref. [14], the desired quantity for an implementation of a DFT+- 1- 
method is a projection 

V c = J2l \L) (L\ of the full Kohn-Sham Green function Gks(w) on a set of localized 
orbitals {|£)}: 

G£sH = V c G KS (cu)V c . (2) 

The subspace C = span({|L)}) is usually termed correlated subspace. It is the subspace 
of orbitals in which many-body correlations play a major role and where corrections to 
the DFT will be considered. In plane- wave based calculations, Gks(w) is available in 
terms of a (truncated) set of Bloch states \K) that are eigenstates of the Kohn-Sham 
Hamiltonian iJxs \K) — e K 

Inserting equation (3) into equation (2) shows that one needs to evaluate projections of 
the type (L\K) in order to access the matrix elements (G^sO-^ll' of the local Green 
function. In most cases the correlated orbitals are d or / orbitals, which are to a good 
approximation localized inside the PAW augmentation spheres. For \L) within these 
spheres and given the PAW decomposition of a Bloch state \K) one obtains 

(L\K) = (L\K) + (L\ - ft)) fr\K) \ , 

which simplifies for a converged partial wave expansion to 

(L\K) = j2(m)mK). 
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The index % of the augmentation functions |0j) includes site /z, angular momentum / 
and m as well as an index v labeling the radial function: i = (/x, l,m, is). In practice, 
the localized orbitals can be also chosen to be angular momentum eigenstates at a given 
site /x, which leads to 

(L\K) = J2(L\<Pu)(Pu\K) (4) 

V 

with v abbreviating i = (yU, /, m, v) where /i, Z, and m are fixed. In the PAW approach, 
the first augmentation function, v = 0, for each channel is usually taken to be an 
atomic eigenfunction (c.f. Ref. [17]). Defining the correlated subspace in terms of 
atomic eigenfunctions leads consequently to \L) = |0„=o) : 

(<P u=0 \K)=J2(<P»=o\<Pu)(Pu\K). (5) 

V 

As higher augmentation functions, v > 0, are in general not orthogonal to the \<j) u =o) 
state, evaluation of equation (4), requires accounting for the overlaps of the form 
(4>u=o\4>iy') ■ This approach has been implemented to the AB-INIT code [20, 21] by 
Amadon [10]. In the present work we have used the Vienna Ab-initio Simulation Package 
(VASP) [18, 19] to implement the same method. It will be referred to as PLO(A) in 
what follows. 

The first order approximation to equation (5) has also been implemented to 
investigate how accurate it is. In this case one only uses the first term of the sum 
in equation (5) 

(<f>„ =0 \K) = (0, =o |0,=o) (Pu=o\K) (6) 

while disregarding all other terms. This approximation will be referred to as PLO(0). In 
addition to the above methods we have implemented another scheme: As in the LDA+U 
scheme implemented in VASP [22, 23] we choose 

\(L\K)\ 2 = J2(K\P„)(^>)(P»>\K). (7) 

U,l/' 

Hence, the absolute value of the projections is in this scheme given by the projection 
onto a subspace of augmentation channels with given angular momentum, (/,m). The 
phase is determined by 

arg((L|^)) = arg^(^|K)j. (8) 

This particular construction will be referred to as PLO(V) in the following. If higher 
augmentation channels are negligible, | {p u= o \K) \ ^> | (p u>0 \K) |, and additionally 
(4>v\<t>v') = o~ij (which is not the case in general) equation (6), equation (5) as well as 
equations (7)- (8) become formally identical. The approach defined in equation (7)- (8) 
differs from the approach in equation(4) in one important point. In the former approach 
a specific radial function |0o) is used for the projection and only overlaps of the type 
(0o 1 are taken into account, whereas in the latter approach the radial dependence 
is averaged out by including general overlaps (0^10^). As constructed, the projections 
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in equation (5) as well as in equations (7) and (8) are not properly normalized for two 
reasons: (1) the Bloch basis is incomplete since only a limited number of Bloch bands is 
included and (2) the PAW augmentation functions are in general not orthonormal. In 
our implementation we orthonormalize the projection matrices by the following Wannier 
type construction: By definition, the localized states \L) are labeled by site and angular 
momentum indices: L = (/j, Z, m). We split the site index // = R + T such that R labels 
the position within the unit cell and T is the Bravais lattice vector of the unit cell in 
which n is located. This allows us to construct the Bloch transform of the localized 
states, 

|£ k > = E e * kT | L T>, (9) 

T 

where k is from the first Brillouin zone and \L T ) = \L) = /, m) with /i — R + T. The 
sum in equation (9) runs over the Bravais lattice. Labeling the Bloch states \K) = |k, n) 
by their crystal momentum, k, and band index, n, we normalize our projection matrices 
Vl n (k) = (Lk|k, n) using the overlap operator 

o w (k) = J2v c L M(v c L/n (k)Y (io) 

n 

in 

^L(k) = £[0(k)]^ 2 P£, n (k). (11) 

v 

These orthonormalized projection matrices are calculated once at the beginning of 
any calculation and can then be used to obtain the local Green function of the correlated 
orbitals from the full Bloch Green function 

G e L M = -P c Ln ^)G^ ch (k,u) (^, n ,(k))*. 

k,rm' 

Similarly the hybridization function, A(u), is available. It is related to the local 
Green function by 

G-\u) = uj + i5-e d - A(w), (12) 

where e d is the static crystal field, equation (12) is a matrix equation with G, A, 
and €d being (dim C) x (dim C) matrices, in general. To separate the hybridization 
from the static DFT crystal field, we numerically evaluate the limit u> — > oo, where 
u - e d . 

In a DFT+DMFT calculation the projection matrices P1 n (k) are used for up- and 
downfolding quantities like the Green function and the self energy in the course of 
the iterative DMFT procedure in exactly the same way as shown for the local Green 
function above. For example, the self energy obtained by an impurity solver for the 
effective impurity model TP LL ,(oo) can be upfolded to the Bloch basis as follows: 

£^r(k,u,) = J2 (Xn(k))* TF LL ,{u) 7* B ,(k). 

LL' 
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Figure 1. DFT(LDA) total and orbital projected density of states for SrVC>3 



Since the self energy in DMFT is a purely local quantity, the index k on S^? cfi- (k, u) 
reflects the momentum dependence brought about by the projection matrices. The 
presented projection scheme allows for the inclusion of both correlated and uncorrelated 
states in the procedure. Therefore, information about the interplay of correlated orbitals 
with their uncorrelated ligands can be obtained. 

With these quantities available we can proceed to perform DFT++ calculations of 
different kinds. The local Green function or the hybridization function can be used as an 
input for impurity calculations, e.g. impurities on surfaces as discussed for the examples 
of Co on Cu(lll) in Ref. [24] as well as Co on graphene in Ref. [25]. The DFT+DMFT 
implementation based on the projectors from equations (7)-(8) has been compared to 
NMTO [13, 26] and other plane wave Wannier function based implementations [27, 16] 
for different cases (e.g for Sr 2 Ru04, LaTiOs, YTiOs, NiO), where the results were in 
good accordance [28, 29]. 

2.2. Many Body Methods, Impurity Solvers 

The goal of the present DFT++ implementation is to interface DFT with a 
subsequent many-body treatment. Since our implementation yields Green functions 
and hybridization functions, we can apply it to any many-body method based on those 
quantities. A widely used DFT++ methodology is the above mentioned DFT+DMFT. 
Here the local downfolded DFT Green function, which contains the full hybridization 
function, serves as an initial input to the impurity solver. We have used here the Hirsch- 
Fye Quantum Monte Carlo [30] (HF-QMC) solver. In all cases an additional input is the 
local Coulomb interaction matrix. The DFT+DMFT Hamiltonian can then be written 
as follows for a two or three band model 




(13) 




+ 2 ^ ] Jmm' (^ c lrn,a C lrn' -a C im-a c im' ,cr 
~l~ C irn,<r C im,-a C irn' ,crCi m ' ■ 

Here ra^ m is the number operator for electrons at site i with spin a in an orbital 
characterized by the quantum number m. The single-electron part given by i^KS contains 
all material specific quantities. 

The terms in the first two lines of equation (13) contain only density- density 
type interactions. In the following lines, however, the so-called spin-flip and pair- 
hopping terms appear in which two electrons in different orbitals flip their spins or 
transfer to other orbitals. If the HF-QMC method is used as an impurity solver, one is 
limited to density- density interactions. The recently developed continuous-time Monte 
Carlo (CT-QMC) [31, 32] impurity solvers can process the full rotationally invariant 
U-matrix. It is important to note that the neglect of the spin-flip and pair-hopping 
terms happens for technical reasons in the first place. It is not always a reasonable 
assumption that the exchange coupling J is small with respect to U and can be left out 
to a good approximation. The neglected terms can be of importance in certain systems 
[33, 24, 34, 29]. If the local part of the Hamiltonian (13) is solved by the QMC method 
continuous hybridization functions can be taken into account. While QMC is formally 
exact, the method provides data with statistical noise on the imaginary axis, which 
makes analytic continuation problematic. 



3. Applications 

3.1. SrV0 3 

3.1.1. DFT The correlated metal SrVOs has become a common testing ground for 
first-principles many-body techniques and has been the subject of multiple theoretical 
and experimental investigations [35, 36, 37, 38, 39, 40, 41, 13, 42, 43, 44, 45, 46, 12]. It 
has a perfectly cubic perovskite structure (space group Pm3m) with a lattice constant of 
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Figure 3. Effective impurity local density of states as obtained by PLO(V) (green), 
PLO(A) (red), and PLO(O) (blue) methods using the V ti g as well as the O 2p bands 
in the projection. 




-8 1 1 1 1 1 

r r x m r 



Figure 4. Difference between the amplitudes of the orthonormalized ti g projectors in 
the PLO(O) and PLO(A) methods for V t2 g as well as the O 2p bands in the projection. 
15 eV-<5„(k) (see equation 14) is shown as the width of the corresponding bands. 

3.84 A [47]. The V ion is surrounded by O ligands in a perfectly octahedral configuration 
leading to a splitting of the d orbitals into t 2g and e g crystal-field states. The density of 
states (figure 1) and band-structure (figure 2) obtained using the VASP code and LDA 
reveals three isolated bands at the Fermi level which are formed by the degenerate t 2g 
states of vanadium. The bandwidth of this block amounts to 2.5 eV. The Van Hove 
singularity leV above the Fermi level in the DOS corresponding to these bands shows 
their predominantly two dimensional character. 

The e g states of vanadium are mostly located above the t2 g bands as conventional 
ligand-field theory suggests [48] . The bands below the V d bands, extending from -8 eV 
to -2 eV are predominantly of oxygen 2p character, but also show hybridization between 
O 2p and the V d states (figure 2). 

The whole series of transition metal oxides SrV03-CaV03-LaV03-YV03 ranging from 
the metal SrV03 to the insulator YVO3 can be classified, following Zaanen et al [77] , as 
Mott-Hubbard systems, since the ligand p bands are clearly separated from the transition 
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Figure 5. Impurity spectral function obtained by DFT+DMFT(QMC) for U = 4 eV 
at the inverse temperature (3 = 20 eV -1 . Only the three V t 2g bands were used in the 
projection. 



metal d states at the Fermi level. Thus, the low energy physics of the material is mainly 
determined by the three V t 2g bands around the Fermi level. This suggests to use 
only these bands for the construction of the effective low energy Hamiltonian for the 
DFT+DMFT calculations. Such a construction, however, will not give any information 
of the behaviour of the O p states, yet this might be of crucial importance for the 
description of the physics of the material. An important example where this is the case 
is NiO, where the charge-transfer behaviour of the system cannot be described by taking 
only the Ni d states into account [49, 50, 51]. 

3.1.2. Tight-Binding discussion of the 3 band case The tight-binding like Hamiltonian 
created via the PLO method from only the V t 2g states is compared to other schemes for 
generating Wannier functions, the maximally localized Wannier functions (MLWF) as 
defined by Marzani and Vanderbilt [15] in the PAW and the FLAPW basis sets and the 
Nth order muffin-tin method (NMTO) [52] in table 1. To compute the hopping matrix 
elements two steps are necessary. First, one obtains the Hamiltonian in the Wannier 
representation by application of the projection matrices to the Bloch Hamiltonian 
H^ n ,(k) = e n (k)5 nn ' in the following way 

tffL'(k) = £^( k )^nn'(k) (^, n ,(k))* 

nn' 

Second, Fourier transformation of the k-dependent Wannier Hamiltonian H^ L , (k) yields 
the on-site energies and hopping integrals 

1 N 

= n £ exp ^ k ' ( R " R '))^(k). 

k 

In above formula L, L' label the orbitals and R, R/ are lattice vectors. 

The nearest neighbour hopping clearly dominates, which shows the short range 
of the bonding in SrVC>3. The three compared methods although very different in 
cost and conception yield virtually identical descriptions of the system. An explicit 
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Figure 6. Momentum resolved impurity spectral function obtained by 
DFT+DMFT(QMC) for U = 4 eV. The three V t 2g bands were used in the pro- 
jection. As a comparison the DFT(LDA) band structure of the V ti g Bloch states is 
shown. 



comparison of the PLO method with MLWFs generated by the WANNIER90 package 
[53] yields quantitative agreement with the present implementation in the present case. 
This underlines the validity of the our implementation. This is another instance of 
the known fact, that the first guess in the procedure of maximal localization in the 
MLWF scheme is usually already a very good one for transition metal oxides [54] . Since 
our implementation is such a first guess procedure without any additional localization 
procedure the good agreement is not surprising. Additionally, the specifics of the system, 
like its high symmetry and the well separated block of t2 9 bands make it a rather 
elementary case. 



3.1.3. Comparison of different PLO methods Differences in the different PLO methods 
can become significant if the oxygen 2p bands are included. These states are essential 
in understanding the physics of many transition metal oxides. Their importance for the 
physics of the SrV03-CaV03-LaV0 3 -YV0 3 series has been pointed out by Mossanek 
et al [55]. The oxygen 2p bands are located below the V ti a block and hybridize 
considerably with them as shown in figure 2. The number of bands to be taken 
into account in the Wannier construction now is 12. The resulting local densities of 
states (LDOS) for the effective three band model are shown in figure 3. While the 
LDOS created by the PLO (A) and PLO(V) methods are virtually indistinguishable, 
the approximation PLO(O) creates a different LDOS. The distribution of the spectral 
weight is different in the PLO(O) method. In fact, the amplitudes of the projectors 
PLO(V) and PLO (A) are close to identical, while the PLO(O) projectors can differ 
significantly from the other two methods, leading to the difference in the resulting 
LDOS. The difference <J n (k) between the projections PLO(O) (PL( k )o) and PLO(A) 
(P£ n (k)^) is shown pictorially in figure 4 for every band n. Here the total difference 
between the amplitudes of the t 2g projectors 

^(k)=E|l^n(k)^|-pL(k)o| 2 (14) 
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has been computed. For better visibility 15 eV-<5 n (k) is shown as the width of the 
corresponding bands. On the same scale the difference between PLO(A) and PLO(V) 
would be by far smaller than the linewidth used for drawing the LDA bands in figure 4. 
The values attained in this case are <5 n (k) < 5.5 • 10 -5 . This figure gives a qualitative 
picture of where the higher order terms of the expansion equation (4) give significant 
contributions. This is the case outside of an energy window extending from —0.5 eV 
to —2.5 eV where no significant difference between PLO(0) and PLO(A) is observed. 
The quite large differences outside this window stem from the fact that the first order 
projector is constructed as the atomic eigenfunction and thus the projected bands from 
it overweight the lower lying bands leading to an underestimation of the spectral weight 
at and above the Fermi level. This leads to a higher occupancy of the effective impurity 
in the PLO(0) approximation as compared to the other methods. The PLO(A) and 
PLO(V) methods yield impurity occupancies of 0.7 electrons per orbital (including 
spin degeneracy), while the PLO(0) yields a filling of 0.85. This stems from the fact, 
that the impurity level is much lower in PLO(0), because the lower lying bands are 
overweighted. In fact, the impurity level (static crystal field) lies at about —1.1 eV 
when using the PLO(0) method, whereas it it is centered at approximately —0.6 eV 
for the other methods. A difference between the PLO(A) and PLO(V) constructions is 
that in the former a specific radial function \<j) v =o) is chosen for the projection, while in 
the latter an averaging over radial dependencies takes place by the inclusion of general 
overlaps {4> v \4> v i) . This did not make any noteworthy difference here or in other systems 
we have considered. 

3.1.4-. DFT+DMFT calculations Let us now turn to the results obtained by 
DFT+DMFT using the projection scheme explained above. We have performed 
calculations using the HF-QMC solver for the material. We obtain momentum resolved 
spectral functions and are thus also able to compare our results with recent ARPES 
studies [42, 56]. For the calculations including only the t 2 g bands we have used the 
on-site interaction U — 4 eV and for the calculations including also the O p bands 
U = 6 eV was used. In all cases J = 0.65 eV was applied. This set of parameters was 
chosen to be able to compare with earlier studies of the material that used the same 
values [14, 10]. Calculations were performed at inverse temperature (3 = 20 eV" 1 using 
200 time slices and 10 6 Monte-Carlo sweeps. Spectral functions were obtained from the 
measured Green function via the maximum entropy method [57]. 

First we examine results obtained by the minimal 3 band model obtained only from 
the three V t 2g bands. As in earlier studies [14, 10, 45, 12] we obtain the local impurity 
spectral function for the V ti g states. The typical three-peak structure is apparent in 
the spectral function in figure 5, with upper and lower Hubbard bands at ~ 3 eV above 
and at ~ 2 eV below the Fermi level, respectively. The data agree very well with other 
studies [14, 10, 45, 12]. We also estimated the mass renormalization 
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Figure 7. Orbitally resolved spectral function obtained by DFT+DMFT(QMC) and 
PLO(0) (left) and PLO(V) (right) for U = 6 eV. The V t 2g and O 2p states were used 
in the projection. 



from the imaginary part of the self energy at the first Matsubara frequency to be 
Z = 0.61, which means, that the mass enhancement is m*/m ~ 1.65 in accordance 
with ARPES estimates of m*/m ~ 1.8 ± 0.2 [42]. We also computed the momentum- 
resolved spectral function from the QMC data 

Ai(k, u) = Im (cu + n - Ei(k) - E^w)) -1 , 

7T 

which is shown in figure 6 together with the three V t 2g DFT(LDA) bands as a 
comparison. The Hubbard bands are clearly recognizable as non-coherent features 
around —2 eV and in the range of 2 eV to 4 eV above the Fermi level. The agreement 
with ARPES energy dispersions is generally quite good, especially the bottom of the 
quasiparticle band is found at —0.6 eV in contrast to the LDA value of —1 eV [42]. The 
overall width of the quasiparticle bands around the Fermi level is reduced to about 1.5 
eV, similarly to earlier reported data [45]. 

Since the inclusion of the oxygen 2p states into the model is important we also show 
results for that situation. The Wannier functions created in this case will be more 
localized than in the t 2 g only model, because they are now composed of a larger number 
of Bloch functions. The effective impurity problem is now constructed encompassing 
also the p bands and thus also the residual d spectral weight contained in them, see 
figure 2. The impurity spectral function will now also have some spectral weight inside 
the oxygen block. The spectral functions obtained via two different PLO methods while 
using QMC as impurity solver within DMFT are shown in figure 7 for PLO(0) and 
PLO(V), respectively. In this case the parameter values U — 6 eV and J = 0.65 eV were 
applied. The impurity clearly shows the spectral weight inside the oxygen block which is 
stronger for the PLO(0) for reasons outlined above. The quasi-particle renormalization 
is Z = 0.62 (from PLO(V), from PLO(0) it is slightly higher at Z = 0.66) in this case, 
which is in agreement with the previous t2g-only estimate. Along with the impurity 
spectral function figure 7 also shows the spectral function decomposed by the respective 
Bloch bands. These can be obtained by applying the inverse of the projection matrices 
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Figure 8. Momentum resolved impurity spectral function obtained by 
DFT+DMFT(QMC) and PLO(V) for U = 6 eV. The V t 2g and the O 2p states 
were used in the projection (12 bands). Again the LDA band structure of the V t2 g 
and O 2p Bloch states is shown for comparison. 

and thus upfolding the impurity back to the Bloch basis. The assignment of a certain 
band to a certain group, say O p is performed using the dominant character of the band, 
as in the DFT analysis. The lowest lying nine bands are thus labeled as oxygen bands, 
the following three bands are labeled V t 2g . The differences between the PLO(O) and 
PLO(V) methods carry over into the DFT+DMFT description. The effective impurity 
shows considerably more spectral weight inside the oxygen block around —5 eV for 
the PLO(O) consistent with the LDA LDOS. While the PLO(A) and PLO(V) yield 
occupancies of 0.66 per orbital one finds 0.81 with PLO(0). The differently distributed 
spectral weight also gives rise to different densities of states at the Fermi level. For 
the PLO(V) case it is ~ 0.62 eV" 1 while it is reduced to ~ 0.56 eV" 1 for PLO(0). 
Furthermore, the quasi-particle peak is considerably reduced in PLO(0) at the expense 
of the enhanced spectral weight at lower energies. The upper and lower Hubbard bands 
also appear shifted in the direction of the Fermi level by some 0.5 eV. These effects 
show, that the PLO(0) approximation does not capture the whole spectral weight of the 
t% g states around the Fermi level correctly. The momentum-resolved spectral function, 
calculated using 12 bands in the Hamiltonian, is shown in figure 8. The t 2g bands 
around the Fermi level are in close agreement with the ones obtained in the t2g-only 
case, showing a considerable renormalization as compared to the LDA bands, as before 
in accord with experimental studies. The upper Hubbard band is visible at 3 eV above 
the Fermi level, while the lower Hubbard band is hidden inside the oxygen bands starting 
at about 2 eV below the Fermi energy. The oxygen bands appear slightly broadened 
and shifted with respect to the t 2g bands, but essentially unchanged as compared to the 
LDA bands. 

The inclusion of oxygen states gives rise to the so-called double counting problem 
of DFT+DMFT, or more generally of any DFT++ method which includes correlated 
and uncorrelated states. 

To avoid a double counting in the energy a term H& c = ^dcJ2ma n m,a, with 
the double counting potential /i^ and orbital index m, has to be subtracted from 
the Hamiltonian, because the local Coulomb interaction is already contained in DFT 
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in an averaged manner. The double-counting amounts to a shift of the chemical 
potential of the correlated orbitals and is an issue, since there is no direct microscopic 
or diagrammatic link between the DFT and the model Hamiltonian part. In the 
DFT+DMFT implementation that we used a constraint on the electronic charge to 
define a double counting correction. The occupancies of the correlated subspace 
computed from the local non-interacting Green function and the interacting impurity 
Green function are required to be identical, which can be stated in the following form 



The double counting potential fid c is iteratively adjusted to fulfill equation (15). The 
most obvious effect of a change in the double counting is a shift of the oxygen states 
versus the V d states. The effect of the double counting is not so dramatic for the case 
at hand (as shown in Ref. [10]), yet it influences quantitative results. This problem can 
be resolved by explicitly calculating the interaction between the d and the p states (U p d) 
or in methods like GW+DMFT [58], where the double counting is explicitly available. 

3.2. Mott transition in Ca 2 ~ x Sr x Ru04 

As a second application we present the Mott transition in Ca 2 _ :r Sr. r Ru04 at small 
Sr contents (x < 0.1). This system is an example of a strongly distorted perovskite 
crystal, where the local symmetry is reduced. We will show that the PLO methodology 
we presented can also be applied to systems with lower symmetry and that we can 
obtain accurate results concerning the coupled structural and electronic metal-insulator 
transition in the system. The system has been studied during recent years, since it shows 
a variety of interesting structural, electronic and magnetic phases [59, 60]. Especially 
the discovery of unconventional superconductivity in the pure Sr compound (x = 2) 
by [78] has caused interest in the correlated electronic structure of these ruthenates. 
Substitution of Sr 2+ with the smaller, but isoelectronic, Ca 2+ results in structural 
distortions of the cubic (IAi/mmm) crystal structure of Sr 2 Ru0 4 . These distortions 
reduce the crystal symmetry to the Pbca symmetry group and transform the system 
into an antiferromagnetic Mott insulator at low temperatures in the doping range (0 < 
x < 0.2) [60]. 

The low energy physics of the material series is governed mainly by the Ru-4d(£ 29 ) 
bands close to the Fermi level [61] that are occupied by 4 electrons. The tetragonal 
symmetry eventually leads to an additional "2D" crystal-field splitting within the t 2g 
bands, breaking their degeneracy [62]. We will refer to this additional splitting as the 
crystal-field (CF) splitting. The crystal structure of Ca 2 Ru0 4 in the low temperature 
S — Pbca phase is shown in figure 9. The structural distortions can be classified as 
rotations of the RuOg octahedra around two different axes and an additional compression 
or elongation of the octahedra along certain of their symmetry described in Ref. 

[60]. The effects on the electronic structure can be understood in terms of hybridization 
of the Ru(4d) orbitals with the neighbouring 0(2p) orbitals as has been investigated in 
Refs. [63, 62]. The coupled structural and electronic metal-insulator transition occurring 




(/3) = itc£5(fl- 



(15) 
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Figure 9. Crystal structure of S — Pbca Ca2Ru04. Only atoms contained in the unit 
cell are shown: Ru (blue), O (red) and Ca (green). The octahedron surrounding every 
Ru atom has O atoms on each of its corners. The Ru layers are separated by rock 
salt-like layers of Calcium and Oxygen. 
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Figure 10. Effective impurity local density of states in the crystal field basis (CFB) of 
the S — Pbca structure of Ca2Ru04 for the effective three band model (blue) compared 
to the full Ru d density of states (red). 



in the material [64] has been a matter of controversy. On increasing the on-site Coulomb 
interaction U Anisimov et al [65] found an orbital selective Mott transition in the ti g 
bands. A different mechanism for the metal-insulator transition was proposed by Liebsch 
and Ishida [79] who investigated the transition with respect to the Coulomb interaction 
U and the CF splitting and found a common metal-insulator transition in all t<i g orbitals. 
In both works local quantum correlations were taken into account correctly. The effects 
of crystal distortions were, however, neglected completely or only incorporated using 
model parameters, like the CF splitting. Since both structural modifications and local 
dynamical correlations govern the physics of the system, it is imperative to include the 
structural distortions in a true ab initio manner. 

3.2.1. DFT For the ruthenate system three Wannier orbitals corresponding to t<i g 
symmetry are used. The construction of those orbitals is complicated by the symmetry- 
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Figure 11. DMFT spectral functions for the metallic structure at (3 = 20 eV _1 (left) 
and for the insulating one (3 — 30 eV _1 (right) showing the evolution towards the 
insulating state for x = 0.0. The Ru(4g?) density of states from DFT(PBE) is shown 
in gray. 



reducing structural distortions in Ca2- : rSr x Ru04. The PAW-projectors are defined 
with respect to an intrinsic (global) cubic coordinate system which coincides with 
the crystallographic axes in the present case. The crystal distortions in the system, 
especially the rotational distortions, locally modify the orientation of the orbitals inside 
the Ca/Sr cube surrounding the RuC>6 octahedra. This leads to the fact, that the global 
symmetry character t 2a vs. e g for the Ru(4d) orbitals is not resolved satisfactory and 
we obtain a mixture of global t 2g and e g at the Fermi level. Since it is known from 
previous calculations [61, 63, 62] that the electronic structure of Ca2- x Sr x Ru04 around 
the Fermi level is governed by bands deriving from the Ru(t2 9 ) states, local rotations 
of the orbitals have to be considered in order to be able to obtain the correct orbitally 
resolved density of states and correct projectors in the PLO framework. In this context 
local means that the orbitals are rotated on each Ru site in the crystal cell. The angular 
parts of the five Ad orbitals are represented by real spherical harmonics K\ m (also called 
cubic harmonics) corresponding to the angular momentum I = 2. These functions are 
derived from the real and imaginary parts of the complex spherical harmonics Yi m (see 
e.g. Ref. [66]). Their transformation properties under rotations are therefore related 
to the well known transformation properties of the Y\ m [67, 68] and were given in e.g. 
Ref. [66]. The Ru(<i) projectors in the energy interval considered (~-2.5 eV to ~+0.5 
eV) should show only few e g character since the cubic CF splitting moves the e g states 
to energies of ~ 1.5eV above the Fermi level. Therefore in a first step the contribution 
of the e g orbitals to the projection was minimized. This was achieved by a rotation of 
the underlying coordinate system. The Wannier Hamiltonian in the rotated basis was 
subsequently diagonalized on-site. The latter yields the so called crystal-field basis. The 
orbitals thus obtained are linear combinations of the rotated t 2g orbitals only. They are 
named in ascending order by their energy as (|1) |2) |3)). In the case of Ca2_a;Sr x Ru04, 
the crystal-field basis retains most of the character of the underlying rotated t 2g orbitals 
in the investigated doping range. For x = 0.0 in the insulating (S — Pbca) structure, 
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e.g., it reads 
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The LDOS obtained from the effective three band model in the CF basis is shown 
in comparison to the total Ru-4o? DOS (the VASP code package and the Perdew- 
Burke-Ernzerhof (PBE) gradient corrected functional [80] were used) in figure(10). 
The agreement concerning the bandwidth and general features of the DOS is good, 
discrepancies stem mostly from the e g states, that have been explicitly projected out 
from the CF basis. The effective three band impurity is correctly occupied by 4 electrons. 

Four different crystal structures (obtained via high resolution powder diffraction by 
[60]) have been considered, to ensure that we can picture the metal- insulator transition 
(MIT) both with respect to the temperature T as well as to Sr content x, i.e. x = 
at 180 K (insulating) and at 400 K (metallic), x — 0.1 at 10 K (insulating) and at 
300 K (metallic) [69]. The temperatures here refer to those at which the structural 
measurements were performed in Refs. [60] and [69]. We will refer to the structures as 
low temperature (LT) and high temperature (HT) for every x. Disorder effects emerging 
from Ca — > Sr substitution have not been taken into account; we always calculated 
Ca 2 Ru0 4 in different crystal structures. 

Table 2. Crystal Field Splittings and root- mean-square hopping integrals (see text) 
for all structures considered. 
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It is insightful to first consider the tight-binding picture obtained from the Wannier 
hamiltonian. The crystal field splitting A = £2,3- £1 for each of the structures considered 
is given in Tab. 2. The values for A 2 i = £ 2 — £1 and A 31 = £3 — £1 are very similar, since 
the orbitals |2) and |3) are very close in energy. For clarity only one value is given in 
Tab. 2. It is apparent that the CF splitting is strongly enhanced for the low temperature 
phases at x = and x = 0.1. This is also reflected in the DFT(PBE) occupancies for 
those structures, namely (n 1: n 2 , n 3 ) = (0.70, 0.65, 0.65) and (0.69, 0.66, 0.65) for the HT 
and (0.79,0.62,0.60) and (0.75,0.63,0.62) for the LT structures at x = and x = 0.1 
respectively. The |1) orbital is considerably lowered in energy, its population increases 
via (|2), |3)) — > |1) charge transfer upon change from the HT to the LT phase. The 
nearest neighbour hopping integrals given in table 2 also show a clear picture. We define 
the root-mean square hopping tl ms = N X N X^(on) ^ as a relative measure of the kinetic 
energyf . The sum runs over all nearest neighbours and all orbitals, the normalization 

% this is in the spirit of Ref. [26] 
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Figure 12. DMFT spectral functions for the metallic structure at (3 = 20 eV -1 (left) 
and for the insulating one (3 — 30 eV -1 (right) showing the evolution towards the 
insulating state for x = 0.1. The Ru(4g?) density of states from DFT(PBE) is shown 
in gray. 



factor is given by the product of the number of orbitals N a and the number of nearest 
neighbours N n . The values for the LT structures at x = and x — 0.1 are significantly 
smaller than those for the HT structures. Our DFT results hence show an increasing 
orbital polarization and electron localization as x decreases and also on the transition 
from HT to LT. This is in accord with previous DFT calculations [63, 62]. 

3.2.2. DFT+DMFT calculations Let us now turn to the DFT+DMFT results. We 
parametrize the on-site Coulomb interaction matrix by U mm = U, U mm i = U — 2J, 
and J m ^m' = J [70], with U = 3.1 eV and J = 0.7 eV as obtained from constrained 
DFT calculations by Pchelkina et al [71] for S^RuO^ Since we employ the HF-QMC 
solver to the resulting 3 orbital model, only density- density type interaction terms are 
included. The neglected terms have important effects on the quantitative description 
of the transition. The results presented therefore show an approximate description of 
the system, see the discussion in Ref. [29]. The computations were performed down 
to a temperature of (3 = 20 (~580 K) for the HT structures and (3 = 30 
(~387 K) for the LT ones. We focus on the paramagnetic metal-insulator transition 
exclusively. The spectral functions, shown in figures 11 and 12 with the Ru(4d) DFT 
DOS as a reference, were obtained from the imaginary time Green functions by analytic 
continuation via the Maximum Entropy Method [57]. A small insulating gap begins to 
form at x = 0.1 and becomes clearly visible at x = in the LT phase. The behaviour 
of orbital |1) is qualitatively different from the orbitals |2), |3). This is also reflected 
in the occupancies within DFT+DMFT. At (3 = 30 ^ for x = 0.1 and x = 0.0 they 
read (0.95,0.52,0.53) and (0.98,0.51,0.51), respectively. Thus, the local correlations 
introduced via the Hubbard U strongly enhance the occupation of the |1) orbital driving 
it towards a band-insulating state. The other orbitals approach half filling and undergo 
a standard Mott-transition. In the high temperature structures the introduction of local 
correlations into the system does not alter the results obtained by DFT qualitatively. 
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The system remains in a correlated metallic state. The DMFT occupation numbers for 
the HT structures at x = and x — 0.1 consequently remain almost exactly at their 
DFT values. 

The calculated spectral functions are in good agreement with photoemission data 
for Ca2RuC>4 [72]. The quasiparticle peak in the high-temperature phase, the Hubbard 
bands at —0.5 eV, and at —3 eV can be identified with corresponding features in the 
experimental spectra. The peak at about —1.5 eV binding energy in the experimental 
data reported in Ref. [72] is most probably a contribution of O 2p states, which are 
not included in our calculations. The small gap of about 0.2 eV is in agreement with 
electrical resistivity[64, 73, 74] and optical conductivity [75] measurements for Ca2RuC>4. 
The orbital order in the S — Pbca phase and its disappearance in the L — Pbca phase 
is also in line with the evolution of x-ray absorption spectra (XAS) across the metal- 
insulator transition for x = and x = 0.09 [72, 76]. The observed mechanism of 
the metal-insulator transition is in accord with the mechanism suggested by [79], but 
not to be reconciled with the orbital selective scenario. An orbital selective Mott 
transition would require successive Mott transitions in orbitals |2), |3) and |1) driven 
by the different widths of the subbands. The bandwidth however, which was assumed 
to be the driving force behind the orbital selective transition is of secondary importance 
in the system. 

4. Summary and Conclusions 

We have presented a general interface between a projector augmented wave based all- 
electron DFT method and many body methods based on Wannier functions obtained 
from a projection on local orbitals. Different schemes to obtain projection matrices 
from PAW calculations have been explored and explicitly compared to other schemes, 
like Nth order muffin-tin orbitals or maximally localized Wannier functions for the cubic 
perovskite SrVC>3. 

We find that care has to be taken to correctly represent all the spectral weight of 
the correlated subspace in the construction of the projected local orbitals. The simplest 
scheme, taking into account only the first PAW projector can be a bad approximation 
and leads to an incomplete description of the system. Only constructions taking 
into account higher PAW projectors are capable of a correct description of correlation 
physics on the LDA as well as on the LDA+- 1- level. Our k-resolved spectral functions 
show dispersive quasiparticle features around the Fermi energy resembling renormalized 
LDA bands as well as incoherent features that can be identified as Hubbard bands. 
Mass renormalization factors are consistent with previous experimental [42] as well as 
theoretical studies [13, 45]. Subsequently we have also studied the distorted perovskite 
system Ca 2 _ :r Sr. I .Ru04 for x < 0.1, where a local symmetry adaption of the PLO 
projection matrices is necessary. In accordance with experimental [64, 73, 74, 72, 76, 75] 
and recent theoretical studies [79, 29], we identify a mechanism for the common 
structural and a metal-insulator phase transition: The structural changes in the system 
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enhance the crystal field splitting leading to a Mott transition in some orbitals and to 
a band insulating behaviour in other orbitals. 

Our DFT++ implementation is very flexible and allows for applications ranging 
from the bulk systems, discussed here, to magnetic nanostructures or isolated correlated 
impurities [24, 25]. 
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